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Abstract. Many 3D IV-body barred models of the Galaxy 
extending beyond the Solar circle are realised by self- 
consistent evolution of various bar unstable axisymmet- 
ric models. The COBE/DIRBE A'-band map, corrected 
for extinction, is used to constrain the location of the ob¬ 
server in these models, assuming a constant mass-to-light 
ratio. The resulting view points in the best matching mod¬ 
els suggest that the inclination angle of the Galactic bar 
relative to the Sun-Galactic centre line is 28° ± 7°. 

Scaling the masses according to the observed radial 
velocity dispersion of M giants in Baade’s Window, several 
models reproduce satisfactorily the kinematics of disc and 
halo stars in the Solar neighbourhood, as well as the disc 
local surface density and scale parameters. These models 
have a face-on bar axis ratio b/a = 0.5 ± 0.1 and a bar 
pattern speed flp = 50 ± 5 km/s/kpc, corresponding to a 
corotation radius of 4.3±0.5 kpc. The HI terminal velocity 
constraints favour models with low disc mass fraction near 
the centre. 

The large microlensing optical depths observed to¬ 
wards the Galactic bulge exclude models with a disc scale 
height h z ^ 250 pc around R = 4 kpc, arguing for a con¬ 
stant thickness Galactic disc. The models also indicate 
that a spiral arm starting at the near end of the bar can 
contribute as much as 0.5 x 10 -6 to the optical depth in 
Baade’s Window. The mass-to-A luminosity ratio of the 
Galactic bulge is probably more than 0.7 (Solar units), 
and if the same ratio applies outside the bar region, then 
the Milky Way should have a maximum disc. 

Key words: Galaxy: structure, kinematics and dynamics 


1. Introduction 

Recent results demonstrate that the Milky Way, as more 
than 2/3 of disc galaxies and as suspected early on by 
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de Vaucouleurs (1964) from a comparison of gas kinemat¬ 
ics towards the Galactic centre and in external galaxies, 
is a barred galaxy with the near side of the bar point¬ 
ing in the first Galactic quadrant. Evidence comes from 
near-IR surface photometry, discrete source counts, gas 
and stellar kinematics and gravitational microlensing (see 
Kuijken 1996 and Gerhard 1996 for reviews). The most 
suggestive data certainly are the COBE/DIRBE near-IR 
maps of the Galactic bulge (Weiland et al. 1994; Dwek et 
al. 1995; Binney et al. 1996 for a non-parametric depro¬ 
jection). Estimates of the angle between the major axis 
of the bar and the Sun-Galactic centre line range from 
10° to 45° (e.g. Stanek et al. 1997). 

Furthermore, the distribution of HII regions, young 
stellar clusters, HI gas, CO clouds and dust betray the 
existence of several Galactic spiral arms (Mihalas & Bin¬ 
ney 1981; Vallee 1995 and reference therein), and external 
spiral galaxies of the same (Sbc) Hubble type as the Milky 
Way have arm-interarm surface mass density ratios of or¬ 
der 2 out to at least 3 disc scale lengths from the centre 
(Rix & Zaritsky 1995). Hence the detailed structure of our 
Galaxy clearly deviates from axisymmetry. 

No self-consistent 3D dynamical barred model of the 
Galaxy including simultaneously stellar, gas and dark 
components and extending beyond the Sun’s Galactocen- 
tric distance has been proposed yet. Existing stellar dy¬ 
namical models are either axisymmetric (Kuijken & Du- 
binski 1995; Durand et al. 1996) and/or restricted to a 
single Galactic component (Zhao 1996 for the COBE-bar; 
Kent 1992 and Kuijken 1995 for an axisymmetric bulge 
model), and gas flow calculations always assume a rigid 
rotating bar potential (Mulder & Liem 1986; Wada 1994; 
Weiner & Sellwood 1996). 

This paper presents the first step of a program aiming 
such a complete model. Many self-consistent pure stellar 
dynamical barred models of the Milky Way are built from 
IV-body evolution of bar unstable axisymmetric models. 
This method, already applied to the Galaxy by Sellwood 
(1985; 1993), naturally takes into account the main dy- 
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namical processes acting in the evolution of real isolated 
galaxies, like those responsible for spontaneous bar forma¬ 
tion and self-sustained spiral structures (e.g. Zhang 1996). 
The IV-body method also proves convenient to cope with a 
dissipative gas component, as will be added to the models 
in the next step and described in a second paper. 

The structure of this paper is as follow: in Sect. 2 
we describe the distinct initial conditions of the various 
components considered in the simulations. In Sect. 3, we 
give some technical informations about the IV-body in¬ 
tegration and present the time evolution of the initial 
models. In Sect. 4 we determine for each evolved model 
the best location of the observer (the Sun) according to 
the COBE/DIRBE near-IR observations of the Galactic 
bulge. In Sect. 5 we fix the velocity scales of the models to 
match the observed stellar velocity dispersion in Baade’s 
Window and discuss some of the resulting absolute model 
properties. 


2. Initial conditions 

Equilibrium or close to equilibrium axisymmetric phase 
space density functions (DF) have been obtained for disc 
galaxies either by applying the strong Jeans theorem 
which states that the DF is an explicit function of at most 
three isolating independent integrals of motion (Kuijken 
& Dubinski 1995; Durand et al. 1996), or by solving the 
hydrodynamical Jeans equations for the velocity moments 
under some arbitrary closure conditions and assuming a 
Gaussian velocity distribution (Hernquist 1993). 

The first method has the advantage to provide ex¬ 
act solutions of the Boltzmann equation, but is in gen¬ 
eral (except for Stackel potentials) limited by the lack of 
an analytical third integral: DFs depending only on the 
two classical integrals, i.e. the total energy and the an¬ 
gular momentum about the symmetry axis (z), always 
have <j‘ 2 Rz = 0, unsuitable for substantially anisotropic 
spheroidal components, and = cr RR , which in the Solar 
neighbourhood is wrong for any evolved stellar population. 
Moreover, generating DFs with imposed mass densities re¬ 
quires the cumbersome inversion of the integral equation 
that relates these two quantities (Kuijken 1995 and refer¬ 
ences therein). 

The second method allows for a larger variety of veloc¬ 
ity ellipsoids, depending on the closure conditions, and is 
also more adapted for specified mass distributions. It was 
therefore retained here, with some modification regarding 
the shape of the velocity distribution. 

The initial models are described in standard astronom¬ 
ical units, assuming that the Galactocentric distance of 
the Sun is R 0 = 8 kpc. These units will serve as reference 
in the evolved models until Sect. 5. 


2.1. Mass distribution 


The initial mass density p in our simulations includes three 
axisymmetric components. 

The first one is an oblate stellar nucleus-spheroid (NS) 
inspired from the model of Sellwood & Sanders (1988): 


Pns(s) = 

M ns ( s/a)P 

(1) 

Aira 3 el 00 1 + ( s/a) p ~ q ' 

where 



s 2 

= R 2 + z 2 /e 2 , 

(2) 

loo 

7r 1" (p + 3)7r 

= -CSC - , 

(3) 


a is a knee radius, e the axis ratio and Mns the integrated 
mass. Setting p = —1.8 and q = —3.3, this mass density 
behaves as s -1 ' 8 for s <C a, in agreement with near-IR 
observations of the Galactic inner kpc (Becklin & Neuge- 
bauer 1968; Matsumoto et al. 1982) if a constant mass- 
to-light ratio is assumed, and as s -3 3 for s a, similar 
to the radial number density decrease of RR Lyrae stars 
(Preston et al. 1991) and of the globular clusters (Zinn 
1985). This component is therefore well suited to repre¬ 
sent the nuclear bulge and the stellar halo. 

The second component is a double exponential stellar 
disc: 


Pd(R,z) 


Mu 

4^pT' exp 


R 

h.R 



( 4 ) 


where Hr, h z and Mu are resp. the scale length, the scale 
height and the integrated mass. This component stands 
here for the Galactic old disc. 

Finally an oblate exponential dark halo (DH) with the 
same axis ratio e as the NS component is added to ensure 
a flat rotation curve at large radii: 


PDH(S) = 8 ^ eXp( ~ s/6) ’ 
where b is the scale length and Muu 


( 5 ) 

the total dark mass. 


2.2. Truncation and flattening 


To minimise the number of particles outside the force grid 
and increase the particle statistics near the centre (at fixed 
number of particles), we softly truncate each component 
multiplying its mass density by: 


f(s) = tanh 


s - R c 
26 


( 6 ) 


where R c is the truncation radius and 6 the width over 
which / falls from 1 to 0. The densities thus vanish on the 
spheroidal surface s = R c . 

In all simulations we set R c = 38 kpc, 5 = 5 kpc and 
e = 0.5. The choice of e is realistic at least for the NS com¬ 
ponent and has the technical advantage over more spher¬ 
ical models to limit the ^-extension of the force grid. 
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2.3. Choice of mass density parameters 

The first row in Table 1 lists the values of the parameters 
a, Mns, h.R , h z , Md, b and Mdh that were adopted for the 
“central” simulation mOO. The choice of a and Mns resp. 
refer to the transition radius where the Galactic near-IR 
emissivity starts to deviate significantly from the inner 
power law, and to a mass-to-light ratio in the nucleus re¬ 
gion of 0.75 Solar units in the /7-band, corresponding to a 
compromise between the values of 0.5 derived by Binney 
et al. (1991) from gas dynamics and 1 in Kent’s (1992) dy¬ 
namical bulge model. The disc scales Hr and h z are based 
on recent determinations (Fux & Martinet 1994; Ruphy et 
al. 1996; Sackett 1997 for a review) and are also consistent 
with near-IR photometry (Kent et al. 1991). 

The last three parameters Md, b and Mdh were ad¬ 
justed to reproduce the Galactic rotation curve beyond 
3 kpc, as illustrated in Fig. 1, while holding the local sur¬ 
face mass density of the disc into a reasonable range. The 
inner peak of the observed rotation curve is an artifact due 
to non-circular gas motion in the bar region and therefore 
does not need to be reproduced in axisymmetric models. 
The outer rotation curve of the resulting initial mOO model 
is fairly flat inside 25 kpc and then starts to decrease be¬ 
cause of the exponential decline of the DH density and the 
mass truncation. 

At R = R 0: the disc and total disc+NS+DH surface 
densities are resp. 48 M 0 /pc 2 (integrated over all z) and 
66 M©/pc 2 within \z\ < 700 pc, and the disc, NS and DH 
volume densities resp. 0.095, 1.3-10 -3 and 0.014 M 0 /pc 3 , 
summing up to give ptot(/? Q , 0) = 0.11 M 0 /pc 3 . Our lo¬ 
cal spheroid volume density exceeds estimates from low 
metallicity and high velocity stars (Bahcall et al. 1983). 
Part of the excess could be attributed to a missing thick 
disc component in the model. 

The initial conditions need not to mimic precisely the 
Milky Way because the bar instability expected during 
time evolution will drastically affect the phase-space dis¬ 
tribution. Since we do not control a priori the issue of a 
simulation, many runs are required to appreciate the ef¬ 
fect of the initial parameters on the final properties of the 
evolved models. For this reason we have realised, in addi¬ 
tion to the simulation mOO, ten further simulations where 
each of the parameters a, Mns, h,R, h z and Md has been 
separately set to a lower and higher value than in mOO. 
The adopted values and nomenclature for each simulation 
are listed in Table 1. The DH parameters were always kept 
at the same values. 

2-4- Isotropic Gaussian velocity distribution 

Preliminary simulations (Fux et al. 1996) with isotropic 
and Gaussian initial velocity distribution for each compo¬ 
nent were performed to check how self-consistent evolu¬ 
tion would rearrange the phase space distribution of such 
simple initial conditions. 



0 10 20 30 

R [kpc] 


Fig. 1 . Rotation curve of the initial mOO model, with the con¬ 
tributions of each component. The squares are tangent point 
circular velocities based on the HI terminal velocities com¬ 
piled by Caldwell & Ostriker (1981), assuming R 0 = 8 kpc and 
Vo = 220 km/s, as indicated by the cross. The almost coincid¬ 
ing full and dotted lines represent resp. the circular velocities 
before and after relaxing the DH component 


Table 1. Mass density parameters of the initial models. Dis¬ 
tances are in kpc and masses in 10 10 Mg. Model mOO is about 
centred in the explored parameter space. Unspecified values in 
the other models are similar to those in mOO. Model mil differ 
from mOO only by the DH kinematics and model ml2 by the 
subsequent symmetry-free time integration. The parameters 
A/ns and Mdh are masses without the truncation by Eq. (6) 

Model a A/ns h u h z Md b Mdh 


mOO 

1.0 

3.0 

2.5 0.25 

4.6 9.1 32.0 

mOl 

2.0 




m02 

0.5 




m03 


1.6 



m04 


5.0 



m05 



2.0 


m06 



4.0 


m07 



0.13 


m08 



0.35 


m09 




4.2 

mlO 




5.0 

mil 



non-rotating DH 

ml2 



no symmetries 


Simultaneous evolution of all mass components first 
led to unacceptable strong expanding rings of overdensity 
in the disc, mainly excited by radial mass oscillations of 
the DH damping out over a dynamical time scale. The DH 
was therefore individually relaxed during 3 Gyr with im¬ 
posed axisymmetry before releasing the other components, 
suppressing indeed most of the subsequent perturbations 
in the disc. The origin of the DH disequilibrium will be 
discussed in Sect. 2.6. 

In pre-relaxed DH simulations, the velocity dispersion 
of the disc becomes spontaneously anisotropic within a 
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few rotation periods, with a planar anisotropy compati¬ 
ble with first order epicycle theory (Binney & Tremaine 
1987) and ratios between the velocity dispersion compo¬ 
nents around R = R 0 similar to those observed in the 
Solar neighbourhood for the old disc (see Fig. 2 in Fux et 
al. 1996 and Table 2). Hence isotropic initial kinematics 
seems to suffice for the disc. 

However, the kinematics of the NS never reaches the 
observed radial velocity dispersion anisotropy of the local 
Galactic halo stars and instead sustains gravity through 
too fast rotation, exceeding half the circular velocity (see 
Fig. 3c and Table 2). The mean rotation velocity can of 
course be arbitrarily reduced by changing the sign of the 
azimuthal velocities of selected particles, but then the ro¬ 
tation velocity dispersion would increase and in turn de¬ 
viate from the local observations. Some radial anisotropy 
should therefore appear already in the initial velocity dis¬ 
tribution of the NS component. 


2.5. Anisotropic velocity dispersion 


To achieve higher degree of radial velocity dispersion 
anisotropy in both the NS and the DH components, we 
have solved the Jeans equations for more general closure 
conditions than just isotropy. 

Following Bacon et al. (1983), we assume that the ve¬ 
locity ellipsoid points everywhere towards the Galactic 
centre, i.e. a 2 g = 0 in spherical coordinates, with a free 
anisotropy parameter (3=1 — o gg /<j 2 r depending only on 
r and of the form: 


f3(r) 



(7) 


where r 0 is a transition radius and /?oo the asymptotic 
anisotropy at large r: (3 oc r for r <C r 0 , and (3 = (3oo for 
r r 0 . We also assume no other streaming motion than 
rotation about the symmetry axis, i.e. u( = vg = 0 , and 
take as boundary conditions a 2 r = a gs = 0 on the mass 
truncation surface. The details of the numerical method, 
which can in fact also handle anisotropy parameter de¬ 
pending on 9 , are presented in Fux (1997). 

The solutions provide v ^ = Tpj? + cr|^, leaving free 
the relative contributions of organised and random veloc¬ 
ity in the </> direction. As a convenient choice, which en¬ 
sures isotropy near the centre and low rotation for r r Q , 
we set: 


a 


2 


r Wgg + 


r 2 vl 


( 8 ) 


where r 0 is the same parameter as in Eq. (7). 

The values of r G and (3 a0 are restricted by the condition 
> 0 everywhere. In particular, on the truncation sur¬ 
face of our initial mass models, this condition imposes an 
upper limit to (3(r ) in the range eR c < r < R c depending 
only on the potential and its first derivatives (see Fig. 2). 
The resulting strongest constraint is (3{eR c ) < 0.65. 


Table 2. Observed properties of the stellar halo (subdwarfs) 
and the old disc in the Solar neighbourhood. E 0 is the total thin 
disc surface density (i.e. including also the young disc). The 
values of the disc velocity dispersion are averages over 2 . Ref¬ 
erences are: (1) Majewski 1993, (2) Wielen 1977, (3) Sackett 
1997, (4) Kuijken & Gilmore 1989 

Reference 


Sub dwarfs: 

(T rr — ORR 

131 ± 6 km/s 

1 


a <j)(f> 

106 ± 6 km/s 

1 


age = a zz 

85 ± 4 km/s 

1 


V 4 , 

37 ± 10 km/s 

1 

Old disc: 

VRR 

48 ± 3 km/s 

2 



29 ± 2 km/s 

2 


(J zz 

25 ± 2 km/s 

2 


h z 

300 ± 25 pc 

3 

Thin disc: 

Eo 

48 ±8 M 0 /pc 2 

4 


2.6. Dark halo disequilibrium and velocity distribution 

Coming back to the DH radial oscillations, single relax¬ 
ation of DHs with Gaussian initial velocity distributions 
but variable radial velocity dispersion anisotropy based 
on the technique outlined above (Sect. 2.5) indicate that 
the oscillations strengthen with the amount of anisotropy. 
From these experiments we inferred that the mass oscil¬ 
lations are in fact a consequence of the Gaussian tails in 
the velocity distribution and the finite DH mass extent: 
a significant fraction of the DH particles, about 10% in 
the isotropic case and increasing with radial anisotropy, 
have velocities which carry them outside the truncation 
surface, unbalancing therefore the inner equilibrium. 

To overcome this problem, the Gaussian velocity dis¬ 
tributions of the extended DH and NS components have 
been replaced by a bounded 3D distribution build upon 
standard Beta distributions. This new “H 3 ”-distribution 
has four parameters, n, A, p and co (one per adjustable 
velocity moment), and is described as a function of its 
reduced variables £, 77 and £ in Appendix A. 

With the substitutions £ = v^/v 0 , p = v r /v e and 
£ = vg/v e , the distribution is bounded in velocity space 
by a spheroidal surface with principal axes aligned with 
the i> 0 , v r and vg axes and of half-length v 0 along and 
v e along v r and vg. The boundary surface is not taken 
as a sphere because particles launched at a given spatial 
position can afford higher velocities in the tangential di¬ 
rection than in the radial direction without escaping the 
system. Indeed, the boundary velocities v 0 and v e can be 
quantified using the integrals of motion. If the velocity of 
a particle is decomposed in an azimuthal component, 
and a meridional component, u m , then its energy writes 
E = «4 + L 2 JR 2 ) + d>, where L z = R ■ is the 

angular momentum about the 2 -axis and $(7?, z) is the 
gravitational potential. Assuming that the mass density 
is bounded by an equipotential surface d>p and imposing 
v m = 0 when a particle reaches this surface (otherwise the 
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particle would cross the surface and escape), the conserva¬ 
tion of energy and angular momentum yield the following 
condition for the confinement of the particle orbit inside 
the system: 

< + Rp R2 R v l < 2 ~ ®( R > z )] » ( 9 ) 

where R v is the maximum cylindrical radius attained by 
the particle. Thus: 



and v 0 > Ve¬ 
in our initial models, since the truncation surface of 
the mass distribution is not an equipotential, we simply 
take the maximum value of the potential on this sur¬ 
face, i.e. d> p = <l>(i? c ,0), and set R p = R c . For the ini¬ 
tial mOO model, the resulting behaviour of v e and v 0 in 
the plane is shown in Fig. 3a. This will not prevent some 
particles to escape but at least considerably reduces the 
DH disequilibrium and ensures well defined values of the 
boundary velocities everywhere inside the truncation sur¬ 
face. 

If k, A < 2, ji, < .3/2 or oj < 1 , the velocity distri¬ 
bution presents unphysical singularities on the boundary 
spheroid which should in principle be avoided. In partic¬ 
ular, assuming > 0 and substituting the requirement 
A > 2 in Eq. (A10) yields: 



which puts a lower limit on the mean rotation velocity, 
and thus further restricts the range of the parameter r 0 . 
For the adopted mass truncation radius R c and flatten¬ 
ing e, this condition unfortunately prevents NS models 
with mean rotation velocity below 100 km/s when ap¬ 
proaching R = R c (see Fig. 3a). Therefore, to reproduce 
the observed rotation of the stellar halo in the Solar neigh¬ 
bourhood with Eq. (8), we decided to violate this condi¬ 
tion for the NS component and accept irregular velocity 
distributions in a minor portion of the R—z plane. Simi¬ 
larly, the requirements p > 3/2 and u> > 1 put constraints 
on the meridional components of the velocity dispersion 
which could be satisfied in all NS models, but not in the 
hot central part of the DHs. 

The DH pre-relaxations have been maintained even 
with the adopted B 3 velocity distribution to minimise any 
persistent DH induced perturbations on the visible com¬ 
ponents. For these components, the Jeans equations were 
solved using the total potential after relaxation, the noise 
in the potential being reduced by averaging the DH mass 
density in time. The DH relaxation does not modify much 
the starting potential (see Fig. 1). 


1 

0.8 

0.6 

Si 

0.4 
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Fig. 2. Full line: anisotropy parameter (5{r) of the NS compo¬ 
nent in all initial models. Dashed line: upper limit for positive 
v ^ on the mass truncation surface (in model mOO) 

2.7. Choice of velocity parameters 

For the NS of simulation mOO, the parameters r 0 and Poo 
were adjusted to reproduce the observed rotation and ve¬ 
locity dispersion of the stellar halo in the Solar neigh¬ 
bourhood given in Table 2. The adopted parameters are 
r 0 = 4.3 kpc and Poo = 0.65, leading to the very close 
to critical P(r) profile displayed in Fig. 2. The resulting 
initial velocity moments of the NS are shown in Fig. 3a. 

The reversal of the of r versus a'C anisotropy at R ss 
13 kpc is consistent with the kinematics of the blue hori¬ 
zontal branch field stars (Sommer-Larsen et al. 1994) and 
other halo stars (Beers & Sommer-Larsen 1995). The az¬ 
imuthal anisotropy at large radii could reflect the fact that 
stars easier escape in the radial direction, and that of r 
must vanish on the boundary of a finite and stationary 
system, whereas < 7 /^ can still be supported by low eccen¬ 
tricity orbits. 

Figure 3b shows the NS kinematics in simulation mOO 
shortly before the formation of the bar. Within R ^ 14 
kpc, the NS looses a part of its radial velocity dispersion 
anisotropy and its rotation velocity raises. These changes 
certainly reflect the extreme nature of our initial condi¬ 
tions, i.e. forcing the radial anisotropy to its maximum. 
Nevertheless, the re-adjusted velocity moments still re¬ 
main much closer to the observations than those of sim¬ 
ulations started with isotropic NS velocity dispersion, as 
illustrated in Fig. 3c. 

For the NS of the other simulations, we simply use 
the same r 0 and Poo values than in simulation mOO. For 
the DHs of simulations mOO-mlO, we set r 0 = b and 
Poo =0.1, compatible with the restriction of Eq. (12). The 
initial conditions for simulation mil are identical with 
those of mOO, except that the DH has the same anisotropy 
P(r) as the NS and no net rotation, i.e. cr^ = v ^ in¬ 
stead of Eq. (8), and hence also present an irregular outer 
^-distribution. As justified by the preliminary Simula- 
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0 10 20 30 


R [kpc] 

Fig. 3. a Initial velocity dispersion and mean rotation veloc¬ 
ity of the NS component at z = 0 in simulation mOO. The 
two upper dotted lines give the meridional (v e ) and azimuthal 
(v 0 ) escape velocities resulting from Eqs. (10) and (11), and 
the lower dotted line the minimum admissible mean rotation 
(tV,min) defined in Eq. (12) for regular velocity distribution. 
The circular velocity (14) is also represented, b Corresponding 
velocity moments at t = 1200 Myr. c Velocity moments at the 
same time in a simulation identical to mOO except that the NS 
component has isotropic initial velocity dispersion 

tions, all discs have isotropic (/? = 0) and Gaussian initial 
velocity distribution, implying a velocity dispersion tightly 
related to their scale heights. 

3. Time evolution 

The simulations mOO-mll were all done imposing 2-fold 
rotational symmetry about the 2 -axis and reflection sym¬ 
metry about the plane 2 = 0, hence reducing the nu¬ 
merical noise of the potential. For comparison, the initial 


mOO model was also integrated without any symmetry, 
providing our last ml2 simulation. Each simulation is 
runned up to t = 5 Gyr, and ouputs of the particle phase- 
space coordinates were realised every 200 Myr, leaving 
325 evolved models to analyse. 

The number of particles is fixed to 10 5 for the NS+disc 
components and 10 5 for the DH. The proportion of par¬ 
ticles in the NS and in the disc is such that the particle 
masses are exactly the same for both components, and 
may therefore change from one simulation to another: e.g. 
30262 NS particles and 69738 disc particles for mOO. 

3.1. N-body code 

The initial models are integrated with the Particle-Mesh 
code described in Pfenniger & Friedli (1993). 

The potential is computed on a cylindrical polar grid 
using the fast Fourier transform technique in the 0 and 2 
dimensions, where the cells are equally spaced. The radial 
spacing of the cells is logarithmic with a linear core to 
avoid an accumulation point at the centre. The short range 
forces are softened by a variable homogeneous ellipsoidal 
kernel with semi-axes set to 1.1 times the respective cell 
dimensions. The adopted grid has 25(i?) x 24(0) x 201 (z) 
cells and extends up to 50 kpc in R and ±20 kpc in z, 
corresponding to a radial resolution of 40 pc at the centre 
and about 1.8 kpc at R 0 . 

The orbits are integrated using the standard leap-frog 
algorithm with a time step of 0.1 Myr, representing a frac¬ 
tion of the crossing time of the central grid mesh in the 
steep NS potential. 

3.2. Model evolution 

Figure 4 illustrates the evolution of our 13 simulations. 
They all lead to the formation of a bar, but not always on 
the same time scale. Models with higher initial values of 
Toomre’s axisynnnetric stability parameter Q in the re¬ 
gion of rising circular velocity need more time to develop 
the bar, in agreement with the work of Athanassoula & 
Sellwood (1986). This is the case for example when com¬ 
paring the simulations m07 and m08: the former starts 
with a much colder disc, i.e a lower value of Q , and indeed 
forms the bar more quickly. 

The size of the bars clearly varies from one simulation 
to the other. However, one must keep in mind here that 
the evolved models can always be separately rescaled to 
look more similar than they do in initial units, hence com¬ 
plicating an objective comparison. All bars finally flatten 
the surrounding radial profile of the disc, the effect being 
particularly marked in simulation m03. 

Some simulations obviously pass through a double bar 
phase, like in m04t4400. The simulations also develop 
transient spiral arms, especially strong during the time of 
bar formation. Such structures would be hard to achieve 
by other means than the iV-body technique. 
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Fig. 4. Face-on surface density evolution in all runs, including 1/3 of all investigated models. The distances are in initial units 
and the contours are spaced with one magnitude interval. Rotation is clockwise. Note the induced asymmetries in run m!2 
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Fig. 5. Pattern speed flp(f) of the bar in simulations mOO, 
m04 and m08 (initial units). The pattern speed in the other 
simulations roughly fall within the gap delimited by the m04 
and m08 solid lines 

Figure 5 shows the bar pattern speed Hp = (t) as a 

function of time in three simulations. The azimuthal angle 
i9 of the bar major axis (in the inertial frame) is calculated 
by diagonalising the I xx , I yy and I xy components of the 
moment of inertia tensor of the NS+disc particles inside 
R < i? e , where R e is the radius extremising the ratio of the 
inferred diagonal components. The noise level in the Qp(t) 
curves depends on the bar strength and on the presence 
of multiple rotating patterns. 

After the first bar rotations, where the bar may 
slow down by up to 15 km/s/kpc/Gyr, Sip (t) is in gen¬ 
eral constant or slowly decreasing with a rate of a few 
km/s/kpc/Gyr. The lower rotation of the DH in simu¬ 
lation mil does not influence the pattern speed of the 
bar: Slp(t) in mil is almost the same as in mOO. Simula¬ 
tion ml2 has a faster Slp(t) decline than mOO, probably 
related to its noisier gravitational potential. 

The presence of a dissipative gas component may alter 
the bar evolution by speeding up its angular velocity and 
causing its dissolution (Friedli & Benz 1993). 

4. Location of the observer 

The Diffuse Infrared Background Experiment (DIRBE) 
onboard the COBE satellite has mapped the full sky in 
ten different bands, ranging from 1.25^m to 240/mr, with a 
resolution of 0?7x 0?7 and a field spacing of approximately 
0?32. The maps at 1.25 (J-band), 2.2 (A^-band), 3.5 and 
4.9/im, dominated by integrated stellar light, clearly show 
an asymmetric boxy/peanut shaped bulge betraying the 
underlying bar (Weiland et al. 1994). 

The A"-band map offers a good compromise between 
maximum stellar emissivity, low extinction by dust and 
small relative contribution of zodiacal light. If a constant 
mass-to-light ratio is assumed throughout the Galaxy, it 


should therefore provide a reliable tracer of the integrated 
mass distribution and thus put severe constraints on dy¬ 
namical models of the Milky Way. 

This section presents a technique to find the view point 
in an IV-body galaxy from where the simulated panorama 
most ressembles a fixed map, and applies the method to 
our models and a deredened version of the COBE A^-bancl 
data. 

4-1. Reduction of the COBE K-band map 

Even if extinction in the infrared is much less than in the 
optical, it is still not negligeable (about 1.5 magnitude in 
AT towards the Galactic centre). 

Thus, following Arendt et al. (1994), the raw AT-bancl 
map has been corrected for foreground dust extinction by 
first building a J—K color excess map under the assump¬ 
tion of constant intrinsic color, taken as the average color 
in the region l < 30° and 10° < b < 15°, and then trans¬ 
forming it into a AT-band optical depth map using the 
redening law of Rieke & Lebofsky (1985). Such a correc¬ 
tion however is not valid for |5| <3° where a significant 
fraction of the integrated light is emitted along the dust 
screen. Hence this low latitude region must be excluded 
in the fitting procedure. 

Finally, the dust subtracted map has been symme- 
triesed in b and converted from its COBE Quadrilater- 
alised Spherical Cube representation to a cartesian grid 
map in Galactic coordinates with A l = A b = 1° pixel 
size, as shown in Fig. 6. This resolution ensures sufficient 
particle statistics per pixel when computing model maps 
without bluring too much the data. 

f.2. Fitting method 

The position of the observer in the models is specified by 
its distance R in initial units from the “Galactic” centre 
and the scale free-angle <p between the line joining himself 
to the centre and the major axis of the bar, with positive 
<p when the bar is leading the observer. We assume that 
the observer lies in the plane of symmetry, i.e. z 0 = 0, and 
that the visible components have a constant mass-to-light 
ratio Tx in the AT-band. Hence model maps of integrated 
light will depend on R, <p and T k- 

The flux per unit solid angle in a given pixel i is esti¬ 
mated from Monte Carlo integration: 

Ni 

F 1 fR,p,T K ) = (T K An)- 1 Y J Jf^^ ( 13 ) 

k=i k e 

where the sum ranges over all NS+disc particles inside 
pixel i, Ni is the total number of such particles, m the 
mass per particle (identical for both visible components), 
Dk the distance of the k th particle relative to the observer, 
e a softening parameter to reduce the statistical noise in¬ 
duced by the closest particles (e 2 = 0.1), and AH the solid 
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Fig. 6. Grid used for the data interpolation and the model flux 
calculation in the COBE-map fitting. The solid lines roughly 
delimit the 100, 200, 300 and 400 pixels which in the models 
have the best particle statistics 


angle sustained by the pixel. If b\ and &2 = b\ + A6 are 
the pixel limits in Galactic latitude, then: 


Afl = A!(sin &2 ^ sin&i). 


(14) 


The three parameters R , <p and Tk are adjusted to 
the COBE/DIRBE AT-band data by minimising the mean 
quadratic relative residual between model and observed 
fluxes corrected for statistical noise, i.e. (see Appendix B): 


n 2 (R , v,r K ) 


(x 2 




2 V pix 

E 


F ° 2 



-1 


(15) 


where 


X 


2 


pix 

E 


(Fi - F?) 


o\2 


(16) 


lVpi x is the number of pixel, u = N pix — 3 the number of 
degree of freedom, F° the observed flux in pixel i, and 
of the variance of F. t . whose best estimate reduces to (see 
Appendix C): 


Only the bulge region |Z| < 30° and 3° < b < 15° 
enclosing the bar and with reliable dust correction, rep¬ 
resenting 720 pixels, is included in the fits. Moreover, to 
exclude pixels with low number of particles and to check 
the consistency of the fitted parameters with respect to 
the selected l — b sub-region, the pixels are sorted by de¬ 
creasing Ni and only the most populated are retained, 
with iVpi x between 100 and 400. The typical l — b regions 
selected this way are indicated in Fig. 6. 

For each model, 7Z 2 is computed on a 2D grid in R and 
;p with resolution A R = 0.1 kpc and A ip = 1°. The pa¬ 
rameter T k does not require an extra dimension, since for 
fixed values of R and ip the optimum T/c, which minimises 
TZ 2 , can be calculated analytically as: 


Y/y min 


(R,<p) 



where the primes indicate that the T k factor in Eqs. (13) 
and (17) is removed. The 7?. 2 -grid is then smoothed by 
averaging 7 Z 2 over the 3x3 nearest grid points and the 
resulting minimum defines our best fit position of the ob¬ 
server. 


4-3. Tests 

The method has been tested with a simple mass model 
consisting in a triaxial multi-normal bulge for the bar and 
an axisymmetric Miyamoto-Nagai (1975) disc: 


Ptest(x,y,z) = 


Mi 


exp 


b 2 M M 2 

47T 


V8 n 3 a x a y a z 

a M R 2 + (% 


1 

~2 

3z h )(a M 


v2\ 1 


H—T + 


■ z b ) 2 


( R 2 + («M 


■ 2b) 2 ) 5 / 2 ^ 


(19) 


Ni / \ 2 

o? = (T jr An ) - 2 g(^ TF ) . (17) 

The data supply much more accurate fluxes than the mod¬ 
els, so that their errors may be neglected. 

The standard % 2 minimisation method was abandoned 
because the y 2 strongly depends on the precision of the 
F^ s. In particular, the y 2 value will systematically in¬ 
crease with the number of particles available to compute 
the model maps (except for perfect models). Thus the 
fitted view points are biased towards regions where the 
model maps are noisier, and the relative quality of differ¬ 
ent models cannot be judged from this indicator alone. 
Furthermore, the fact that the er 2 ’s are only approxima¬ 
tions of the true variances causes an overestimation of 
the x 2 increasing with the uncertainties on the variances, 
which biases the solutions in the opposite way. This bias 
persists even in 7 Z 2 adjustments since TZ 2 depends explic¬ 
itly on x 2 . To reduce its effect, the of’s have been averaged 
over the neighbouring pixels. 


where R 2 = x 2 + y 2 , z 2 = z 2 + 6^, and with a x = 1.5, 
<j y = 0.8, <j z = 0.4, om = 5 and 6 m = 0.3. The Miyamoto- 
Nagai component is truncated at R = 15 and its mass 
within this radius is 3Mi. 

COBE-like maps were computed for several input ob¬ 
serving points (i?in,^i n ) in the test model and 100 dis¬ 
tinct discrete realisations of this model with 10 5 particles, 
as in the NS+disc components of the dynamical models, 
were fitted to these maps by 7T 2 minimisation, yielding 
the output observing points (i? 0 ut, 7?out)- The statistics of 
the results is shown in Fig. 7. 

For ipi n ^ 40° and Ri n ^ 10, the systematic and sta¬ 
tistical errors on <p out are at most +2° and 5° (except at 
small R[ n for 7V p i x = 100). For N p i x > 200, the statistical 
error on 7? out is less than 4%, and R ou t is underestimated 
with a systematic error increasing with N pix , which is in¬ 
timately connected with the y 2 bias already mentioned 
in Sect. 4.2: the additional pixels in the low flux region 
have the most uncertain a 2 ’s and thus artificially enhance 
the value of 1Z 2 , as indicated in the top plots of Fig. 7. 
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Fig. 7. Test results of the TZ 2 minimisation technique. Ri n and 
i^in are the correct input coordinates of the observer and R ou t 
and y> out the corresponding adjusted output parameters. The 
left and right columns are resp. for variable R ln and fixed y>i n , 
and for fixed R ln and variable y>i n . From top to bottom, based 
on 100 experiments for each input position: the averages and 
standard deviations of 1Z 2 at the minimum, Rout /Rin and ip ou t. 
The full, long-dashed, short-dashed and dotted lines stand resp. 
for iVp ix = 100, 200, 300 and 400 

4-4- Application 

The 1Z 2 method has been applied to all dynamical models, 
noting ip 0 and R a the best fit parameters of the observer 
and 7Z/ 2 uln the value of TZ 2 at the minimum. An example 
of 1Z 2 topography is given in Fig. 8. 

As a compromise between maximum exploitation of 
the data and minimum systematic errors on the fitted view 
point parameters, we have privileged the case N p i x = 300 
and sorted the models by increasing 7^^ nin (300). This 
model sequence was then truncated to include at least two 
models of each simulation, resulting in a sample of 68 mod¬ 
els with 7^^ in (300) < 0.756% (hereafter our “A” sample). 

Figure 9 shows the distribution of the adjusted angles 
ifio for the models belonging to this sample as a function of 



vU 


Fig. 8. 1Z 2 solution for the location of the observer in model 
m08t3200, for N p i x = 300. The 1Z 2 contours are spaced by 
0.5% and the cross indicates the position of the minimum. The 
parameter T k has been substituted using Eq. (18) 

iVpi x . Taking into account the uncertainties of the individ¬ 
ual angles, due to the moment of inertia method for deter¬ 
mining the absolute angular position of the bar major axis 
and to the intrinsic standard deviation of the 77 2 method, 
the 40 best models support an angle of ip Q = 28° ± 7° 
for the Galactic bar, with a possible overestimation by 2° 
according to the tests. 

This result is consistent with the value of 20° — 3CP 
found recently by the OGLE team from the color magni¬ 
tude diagram of red clump giants in several fields across 
the Galactic bulge (Stanek et al. 1997) and with the upper 
limit of 30° set by the MACHO microlensing constraints 
(Gyuk 1996), and confirms the privileged 25° obtained by 
the deprojection of the properly dust subtracted L-bancl 
COBE map (Binney et al. 1996; Bissantz et al. 1996). 

The best models regarding the COBE constraints come 
from simulation m08, started with the thickest and hottest 
disc. Indeed five of the six lowest residual models are from 
this simulation, with 7?.^ in (300) « 0.3—0.4%. For compar¬ 
ison, the initial mOO model has 7 Z 2 > 2.4% everywhere. 
Nevertheless, we note that model m06t4600 gives the best 
COBE-data match within the entire |/| < 90° region. This 
model has a low disc scale height, close to the 204 pc de¬ 
rived by Kent et al. (1991) from the Spacelab IR, Telescope 
data with the assuption of constant h z . 

Table 3 lists the results for a selection of A sample 
models (hereafter our “B” sample) which, in addition to 
the COBE constraints, also best reproduce observations 
in the Solar neighbourhood (see Sect. 5.1). The simula¬ 
tion mOO provide no acceptable model to enter this sub¬ 
sample. 

The distance scale Sr of the models now directly fol¬ 
lows by setting the best-fit Galactocentric distance of the 
observer for N pix = 300 to R 0 = 8 kpc. Figure 10 presents 
the NS+disc face-on configurations of some B sample 
models after such rescaling and the confrontation of these 
models to the COBE data. The B sample models have a 
face-on surface density axis ratio b/a = 0.5 ±0.1 for the 
contours crossing the bar major axis at R = 2 — 3 kpc. 
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Fig. 9. Distribution of the best-fit position angle < p a of the bar in the dynamical models. Dashed line: including all 68 models 
of the A sample; full line: including only the 40 lowest IZmin models of this sample 


Table 3. Location of the observer in the final selection of models (B sample), as constrained by the COBE/DIRBE deredened 
A'-band map. The models are sorted by increasing 1Zmin(N P ix = 300). The last two columns give the adopted distance (Sr) and 
velocity (Sv) scales based on the N plx = 300 solutions 


Model 

JVpix = 
Ro <po[°; 

200 

] RLnl%] 

Ro 

-^Vpix — 
Vo [°] 

300 

R-Ln [%] 

Ro 

AC p i x = 
Vo [°] 

400 

R-Ln [%] 

Sr 

Sv 

m08t4000 

9.3 

31 

0.228 

9.5 

33 

0.342 

9.4 

38 

0.359 

0.842 

0.900 

ml2t2000 

8.2 

23 

0.319 

8.2 

23 

0.409 

8.2 

23 

0.337 

0.976 

0.939 

m06t4600 

9.6 

26 

0.388 

9.6 

26 

0.432 

9.6 

26 

0.619 

0.833 

1.120 

mllt2000 

8.1 

27 

0.365 

9.0 

34 

0.467 

8.2 

28 

0.526 

0.889 

0.964 

m08t3200 

8.6 

24 

0.424 

9.0 

25 

0.472 

8.7 

25 

0.515 

0.889 

0.977 

m09tl600 

8.5 

26 

0.341 

8.5 

26 

0.586 

8.5 

26 

0.828 

0.941 

0.991 

m02t2000 

8.1 

21 

0.547 

8.1 

21 

0.628 

8.0 

23 

0.672 

0.988 

0.946 

m03tl000 

7.5 

22 

0.599 

9.5 

29 

0.641 

8.0 

21 

0.728 

0.842 

1.056 

ml2tl600 

7.3 

20 

0.557 

7.3 

19 

0.677 

7.5 

19 

0.780 

1.096 

1.001 

ml0tl400 

7.J, 

19 

0.580 

8.1 

22 

0.709 

7.4 

15 

0.730 

0.988 

0.916 

m04t3000 

7.9 

18 

0.615 

7.9 

32 

0.738 

7.7 

25 

0.784 

1.013 

0.852 

m06t4800 

8.0 

26 

0.675 

8.0 

25 

0.747 

9.3 

9 

0.650 

1.000 

1.137 


Table 4. Several absolute properties of the rescaled B sample models. X\oc'- mean square residual between model and observed 
local properties; a D ,a s \ local disc and NS velocity dispersions [km/s]; v$ s : local NS rotation velocity [km/sj; h.R,h z : disc scale 
length [kpc] and local scale height [pc]; Ef,pf: local disc surface density [Mg/pc 2 ] and NS volume density [10~ 3 Mg/pc 3 ]; 
Vo- local circular velocity [km/s]; pattern speed of the bar [km/s/kpc] and corotation radius [kpc]; p 0 : bar inclination 

angle [°]; T k- mass-to-K band luminosity ratio [Mq /Lk,q\, to: microlensing optical depth towards Baade’s Window [ICE 6 ]. 
Values in brackets include the DH lenses. The boldfaced models reasonably agree with most of the considered observations 


Model 

Xloc 

D D 
ORR^^ 

ag 

s 

°RR 

s 

erf. 

v^ s 

6r 

h z 

E? 

Po 

Vo 

flp Rl 

9^0 

T K 

TO 

m08t4000 

2.98 

43 

28 

21 

106 

117 

78 

51 

2.7 

363 

47 

0.9 

195 

36 

5.5 

33 

0.60 

1.53 (1.75) 

ml2t2000 

1.91 

44 

27 

19 

118 

107 

89 

63 

4.9 

339 

49 

1.2 

200 

44 

4.5 

23 

0.68 

1.70 (2.12) 

m06t4600 

5.59 

36 

26 

17 

135 

129 

95 

41 

4.0 

232 

58 

1.3 

230 

54 

3.9 

26 

0.59 

1.45 (1.61) 

mllt2000 

1.88 

53 

26 

20 

117 

116 

83 

61 

2.7 

287 

40 

1.2 

210 

52 

4.0 

34 

0.57 

1.38 (1.48) 

m08t3200 

3.88 

42 

25 

22 

111 

124 

81 

58 

2.7 

379 

44 

1.1 

214 

44 

4.8 

25 

0.72 

1.80 (1.97) 

m09tl600 

3.61 

38 

24 

18 

124 

121 

84 

61 

4.4 

259 

50 

1.2 

212 

55 

3.9 

26 

0.73 

2.19 (2.33) 

m02t2000 

1.64 

51 

29 

20 

115 

117 

94 

60 

3.1 

290 

43 

1.1 

210 

56 

3.9 

21 

0.73 

1.97 (2.07) 

m03tl000 

2.47 

41 

30 

18 

129 

125 

93 

44 

4.9 

271 

45 

0.6 

216 

58 

3.9 

29 

0.53 

1.54 (1.89) 

ml2tl600 

3.43 

50 

30 

25 

118 

115 

96 

76 

3.9 

328 

70 

1.7 

217 

48 

4.6 

19 

0.96 

3.00 (3.13) 

ml0tl400 

1.82 

54 

28 

18 

119 

111 

91 

64 

5.1 

295 

47 

1.2 

198 

44 

4.9 

22 

0.70 

1.97 (2.18) 

m04t3000 

3.40 

36 

24 

19 

111 

106 

84 

50 

3.0 

285 

39 

1.5 

200 

48 

4.3 

32 

0.72 

2.03 (2.11) 

m06t4800 

2.39 

47 

27 

23 

133 

133 104 

50 

4.6 

278 

58 

1.7 

229 

45 

4.6 

25 

0.79 

2.28 (2.46) 
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Fig. 10. Comparison of a selection of rescaled models with observations in the bar region. Left: face-on surface density of the 
visible components; The symbol © indicates the location of the observer and the crosses the positions of the Lagrangian points 
Li and L 2 . Rotation is clockwise. Middle: COBE/DIRBE deredened K-band contours (full lines) and corresponding model 
contours (dashed lines). The b scale is dilated relative to l, hence amplifying the deviations. The region below the horizontal line 
was excluded in the adjustments because of unreliable extinction correction. The spacing between the contours is 0.5 magnitude 
in both frames. Right: longitude-velocity diagram of the Galactic HI averaged over |6| < 1.25° and, superposed, the traces of 
non-self intersecting model m-orbits (dashed lines); the innermost trace represents the cusped xi orbit. The model ml2t2000 
has been symmetrised in all frames. The HI data were kindly provided by H. Liszt and refer to Burton & Liszt (1978) 
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5. Model properties and discussion 

Given the distance scale, there remains only one scale 
left in the models since the constraint of virial equi¬ 
librium links together the distance, velocity and mass 
scales. We choose here to fix the velocity scale by re¬ 
quiring that the line-of-sight velocity dispersion of the 
combined NS+disc components towards Baade’s Window 
(, l,b ) = (0?9, — 3?9) equals 113 ± 6 km/s, as observed for 
late-M giants (Sharpies et al. 1990), leaving about 10% 
uncertainty on the absolute mass scale. The resulting ve¬ 
locity scales of the B sample models are listed in Table 3. 

The complete scaling of the models allows now to de¬ 
rive several of their absolute properties, which are pre¬ 
sented and discussed in this section. 


5.1. Local properties 

The main local properties of the models, i.e. the proper¬ 
ties measured at the location of the observer, are given 
in Table 4. These include the velocity dispersion of the 
disc {a° Rl cr^, ofj and of the NS (ag R , cr^ofj, the NS 
mean rotation velocity (v^ s ), the disc scale length (ha) 
and scale height (h z ), the disc surface density (E^), the 
NS volume density (pf) and the circular velocity (To). 

All quantities, except Hr and V c , are based on the 
particles inside a sphere (for the NS) or a vertical cylinder 
(for the disc) centred on the Sun’s position, and hence are 
not azimuthal averages. The disc scale length is an effec¬ 
tive exponential scale length derived from all disc particles 
within the annulus 0.7 < R/R a < 1.2. The circular veloc¬ 
ity rests on the azimuthally symmetrised potential. 





14 


R. Fux: 3D self-consistent IV-body barred models of the Milky Way 


To quantify the ability of the models to reproduce lo¬ 
cal observations, Table 4 also gives the mean y 2 residual 
between the local model properties, from cr RR to and 
excluding h R , and the corresponding observed values men¬ 
tioned in Table 2, weighting the contribution of the three 
disc velocity dispersion components by 1/3 and the four 
NS velocity moments by 1/4. Models from simulation m07 
have very high Xj 2 oc due to their small disc scale height, 
as well as unrealistic mass-to-light ratios and too low mi- 
crolensing optical depths in Baade’s Window (see Sect. 5.3 
and 5.4), and were therefore rejected. 

Many models have disc kinematics consistent with ob¬ 
servations, and in particular correct ratios between the 
components of the velocity dispersion. There appears a 
clear correlation between the vertical velocity dispersion 
of the disc and its scale height, as expected from self- 
gravitating isothermal sheets, although not directly de¬ 
pending on the disc surface density. A proper analysis of 
the disc vertical dynamics would of course require simula¬ 
tions with higher ^-resolution. The disc velocity dispersion 
in model m06t4600 is lower than in m06t4800 because the 
observer is located farther from the centre (in initial units) 
and the velocity dispersion decreases outwards. 

The velocity dispersion of the NS components is in 
general similar to that of the local Galactic halo stars, ex¬ 
cept that its radial anisotropy is insufficient. This differ¬ 
ence is probably linked to the initial model truncation at 
R c = 38 kpc, which limits the amount of radial anisotropy 
inside the models, and the flatness (e = 0.5) of the outer 
mass distribution, which requires substantial support by 
azimuthal motion, either in systematic or random form. It 
is however worthwhile to mention that the local subdwarf 
kinematics in Table 2 may be biased by selection effects. 
For example, kinematically selected samples have much 
larger radial anisotropy than samples selected by metal- 
licity criterions, the velocity dispersion in the latter being 
even consistent with meridional isotropy (Norris 1987). 
Moreover, the kinematics of the globular cluster system 
shows no significant deviation from isotropy. 

All models favour V Q < 220 km/s except those from 
simulation m06. These models have in fact increasing ro¬ 
tation curve from R 0 /2 to i? 0 , so that our velocity scaling 
mainly based on the inner region leads to larger circular 
velocities at R 0 . 

5.2. Corotation and terminal velocities from x\ orbits 

An approximate bar pattern speed Op is derived for each 
model as the least square slope of the bar position angle 
$(t) relative to a fixed axis, within a time interval of a few 
bar rotation and centred on the current time. The related 
corotation radius Rl then follows from the Lagrangian 
points L\ and L 2 , i.e. the saddle points of the effective 
potential 4> e ff = $(i?, z) — |i? 2 Op. Figure 11 plots the 
sequence generated by our models in the Op —i?L plane. 
The B sample values are quantified in Table 2 and the left 
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Fig. 11. Corotation radius Rl, derived as the distance of the 
Lagrangian points L\ and L 2 to the centre, versus the pattern 
speed of the bar Op for the A sample models. Each symbol 
refers to models of the same simulation 


frames in Fig. 10 highlight the position of the Lagrangian 
points. The B sample models have Op = 50 ± 5 km/s/kpc 
and 7 ?l = 4.3 ± 0.5 kpc. 

These results conflict with the dynamical properties 
Op = 63 km/s/kpc and Rl = 2.4 ± 0.5 kpc reported by 
Binney et al. (1991) and are closer to Op = 55 km/s/kpc 
and 7 ?l = 3.6 kpc inferred by Weiner & Sellwood (1996) 
from SPH gas flow modelling in a rigid potential. The 
former Op— R^ pair moreover falls well below the relation 
displayed in Fig. 11. 

Binney et al. (1996) suggest that the peak emission 
in the COBE near-IR maps at b = 0 and near l ~ —22°, 
corresponding to R « 3 kpc, is due to stars trapped by the 
Lagrangian point L 4 /L 5 . We think that this peak could 
also arise from star formation in a gaseous pseudo ring, 
possibly associated with the “expanding 3 kpc arm”, seen 
tangentially and located close to the inner 4:1 resonance, 
as the inner rings in many external barred galaxies (Buta 
1996 and references therein). According to this picture, 
corotation would lie significantly beyond 3 kpc. 

Moreover, Zhao (1996) has constructed a model for 
the Galactic bar using Schwarzschild’s technique and im¬ 
posing flp = 60 km/s/kpc. An A-body evolution of his 
model shows that this input pattern speed decreases by 
about 20% after one rotation period (see his Fig. 9), thus 
becoming consistent with our results. 

Binney et al. (1991) have argued that prograde gas in 
a barred potential, due to the pressure and viscous forces, 
settles onto the stable and closed orbits of the X\ family, 
switching progressively to ever lower energy orbits as long 
as these develop no loops. Depending on the bar poten¬ 
tial, the gas may finally reach a last non-self intersecting, 
“cusped” Xi orbit beyond which shocks transform most 
of the atomic gas into molecular gas, and force the gas 
to plunge to the more viable orbits of the X 2 family (if 
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an inner Lindblad resonance exists). Thus, atomic gas is 
expected to move along non-self intersecting x\ orbits, 
providing the key for a model comparison with Galactic 
HI observations. In particular, the envelop defined by the 
traces of such orbits in the longitude-velocity l—V diagram 
should coincide with the observed HI terminal velocities. 

Orbits of the Xi family have been computed in almost 
all our dynamical models, using the instantaneous frozen 
potential. Figure 10 shows the comparison of several 
model x\ traces versus HI contours in the l — V diagram. 
Some models reproduce fairly well the HI terminal ve¬ 
locity envelop. The best cases certainly are m04t3000 and 
m06t4600, which share the common property to arise from 
simulations with lower disc mass fraction in the bar re¬ 
gion: their disc to total mass ratio within s < 3 kpc is 
less than 0.45, where s is the variable defined in Eq. (2). 
Further cases from the B sample are ml2t2000, m08t3200, 
m09tl600, m02t2000 and m06t4800. All these models, ex¬ 
cept m06t4800, occur shortly after the formation of the 
bar, when redistribution of angular momentum has not 
completely disturbed the initial exponential radial profile 
of the disc. In general, the presence of the bar tends to 
steepen the azimuthally averaged inner radial profile of 
the disc, increasing the central disc surface density. 

Most of the models also have x\ envelopes exceed¬ 
ing the observed terminal velocities, indicating that there 
could be too much mass near the centre. Such an excess 
could of course be reduced by increasing the angle tp 0 , but 
then the peaks of maximum and minimum velocity traced 
by the cusped x\ orbit in the l — V diagram would also 
be shifted towards higher /|. Reducing the velocity scale 
in general render the local disc kinematics less consistent 
with observations. 

Hence, if the Galactic gas really moves on non self- 
intersecting closed orbits, then the HI observations suggest 
that our dynamical models could have too much mass in 
the disc near the centre. 

5.3. Mass-to-K luminosity ratio 

In addition to the position of the observer, the COBE- 
fits in Sect. 4 also yield the A'-band mass-to-light ratios 
T k of the models, without the contribution of the DH 
component. 

Taking Tk,q = 2.69 • 10~ 12 MqW _ 1 Hz (Wamsteker 
1981), the rescaled values for the B sample models range 
from T k = 0.53 to 0.79 in Solar units, except for model 
ml2tl600 (see Table 2). The total NS+disc+DH mass of 
the B sample models within the spheroid s < 3 kpc is 
(2.6 ± 0.15) x 10 10 Mq, with about 8% contribution from 
the DH component. This gives an average percentage by 
which Yk is underestimated. 

Modelling the gas dynamics of the barred galaxy M100 
(NGC4321), which has a Hubble type similar to the Milky 
Way, Knapen et al. (1995) have inferred T k > 0.7 out¬ 
side its nuclear ring. Moreover, for stellar populations with 
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Fig. 12. Confrontation of a model with small h z /R 0 to the 
COBE data. The caption is similar to Fig. 10 

near-IR emission dominated by late K and M giants, as for 
the Galactic bulge (Arendt et al. 1994), similar mass-to- 
light ratios are expected (Worthey 1994). Such lower limit 
for T k rules out some models, like those of simulation m07 
which all have T k < 0.4. The substantial microlensing 
optical depths towards the Galactic bulge also argue for a 
large value of T k, as discussed below in Sect. 5.4. 

A noticeable difference between the COBE A'-band 
map and the model maps in Fig. 10 is that the model con¬ 
tours are steeper in the low latitude and \l\ <; 15° region, 
i.e. where the disc becomes dominant. Even if our correc¬ 
tion for extinction fails at b ^ 3°, the near-IR contours 
in this region still remain very flat after more elaborated 
dust subtraction of the COBE data (Spergel et al. 1997) 
and hence the difference is probably real. Several reasons 
may lead to this departure. 

First, the Galactic disc scale height may be lower by 
a factor ~ 2 in the inner Galaxy than in the Solar neigh¬ 
bourhood, as deduced by Kent et al. (1991) from the 
Spacelab IR Telescope data and by Binney et al. (1996) 
from the deprojection of the COBE A-band map. At con¬ 
stant surface density, discs with smaller scale height are 
more concentrated towards the plane and should there¬ 
fore contribute more to the low latitude near-IR emission. 
Our models, started with radially constant disc thickness, 
do not present such large disc scale height gradient out¬ 
side the bar region. But models from simulation m07 (see 
Fig. 12), with the thinnest disc, indeed have flatter l — b 
contours than models from m08. The variable scale height 
alternative however is not supported by observations in ex¬ 
ternal late-type spiral galaxies (de Grijs & van der Kruit 
1996 and reference therein). 

Second, discs with more foreground mass between the 
bulge and the observer will also enhance the low lati¬ 
tude integrated light. This is the case for the m06 models 
(Fig. 10). At fixed total disc mass, because of the higher 
initial disc scale length, these models have higher disc sur- 
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Fig. 13. a Dependence of the microlensing optical depths to (without DH contribution) towards Baade’s Window on the 
corotation radius Ro,, which is a reasonable indicator of the radial extension of the bar, and on the bar inclination angle <p 0 , 
for the A sample models. The encircled points indicate models with prominent spiral structures, b Relation between to and the 
mass-to-light ratio T k for the same models. The symbols are as in Fig. 11 


face density in the outer region, including the first few kpc 
below R 0 . Increasing simply the disc mass is less efficient 
because the inner disc, mixing with the NS component 
during the bar instability, will also contribute more to the 
bulge emission. 

Finally, the assumption Tx = const, may not hold 
close to the Galactic plane because of near-IR emission 
from interstellar matter, or recent star formation. In the 
latter case, red supergiants, too young for having diffused 
as far above the plane as the old disc population, could 
significantly contribute to the //-emission and reduce T k 
at low Galactic latitude. However our COBE-adjustments 
should not be affected by such a T k gradient since the 
region |6| < 3° was excluded. Alternatively, the IR mass- 
to-light ratio could also be lower in the Galactic disc than 
in the bulge, as found by Verdes-Montenegro et al. (1995) 
for NGC7217 in the /-band. 

In the next section, we will argue for the massive disc 
possibility and against a variable disc scale height. 

5-4■ Bulge microlensing 

The microlensing optical depths r towards the galactic 
bulge provide a serious direct probe of the mass distribu¬ 
tion inside the Solar circle. 

The first results obtained by the OGLE and MACHO 
experiments favour surprisingly large optical depths: the 
reported values are at least (3.3 ± 1.2) x 10 -6 for 9 bulge 
stars essentially located in Baade’s Window (Udalski et al. 
1994), and (3.91)}’®) x 10 -6 for 13 clump giants with mean 
Galactic coordinates l = 2?55 and b = —3?64 (Alcock 
et al. 1997), whereas axisymmetric models predict values 
less than 10 -6 (Evans 1994 and references therein). The 


contribution of bulge lenses may however significantly in¬ 
crease the model predictions if the bulge is an elongated 
bar seen nearly end-on (Evans 1994, Kiraga & Paczynsky 
1994, Zhao et al. 1996, Zhao & Mao 1996). 

We have computed the model microlensing optical 
depths towards Baade’s Window, averaged over all sources 
along the line-of-sight, from the following Monte Carlo 
version of Eq. (5) given by Kiraga & Paczynski (1994): 


T 0S = 


47rG 
Aflc 2 






D s 


( 20 ) 


where the outer sums involve all source particles within 
a solid angle AH of the selected direction and the inner 
sum all lens particles in the same solid angle between the 
specified source and the observer, the symbols D s and D c j 
referring resp. to the source and lens distances relative to 
the observer, and m<i to the mass of the lens particles. G is 
the gravitational constant, c the speed of light and (3 a the 
parameter introduced by Kiraga & Paczynski (1994) to 
describe the detection probability of the sources: (3 S = 0 
if the sources are detectable whatever their distance, like 
clump giant stars, and l3 s « — 1 for main sequence stars 
(Bissantz et al. 1996 and references therein). 

Optical depths are calculated with and without includ¬ 
ing DH lenses. Table 4 gives some results for (3 S = 0. The 
optical depths for other values of /3 S are almost propor¬ 
tional to To: in particular, we find r_i « |r 0 , both with 
or without DH lenses. The DH contributes roughly 9% 
to the total optical depth. In the mean direction of the 
Alcock et al. (1997) fields, 7o is on the average 15% higher 
than in Baade’s Window, in agreement with the above 
reported observations. 
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h z (R„/2) [pc] 

Fig. 14. Microlensing optical depths to (without DH) towards 
Baade’s Window versus disc scale height half way between the 
observer and the Galactic centre for the A sample models. The 
upper frame shows that the correlation is not induced by a 
variable disc surface density. The symbols are as in Fig. 11 

Figure 13a shows the optical depths of the A sample 
models towards Baade’s Window as a function of the size 
of the bar and its inclination angle relative to the observer. 
The optical depth strongly depends on (p 0 , increasing from 
tq ~ 10~ 6 for ip 0 = 40° to To ~ 2.5 x 10~ 6 for ip 0 = 15°, 
and is roughly proportional to Rl. The dispersion in each 
plot partly reflects the dependence of the optical depth on 
both ip 0 and R i/. at fixed Rl, the larger values of To come 
from models with more end-on bars. 

Clearly, the optical depths depends on the mass scal¬ 
ing and thus on the adopted mass-to-light ratio T k, as 
depicted in Fig. 13b. Within the various investigated mod¬ 
els, this relation does not seem to depend much on a third 
parameter. According to it, the ler lower limit of the ob¬ 
served optical depths, r ^ 2.1 x 10 -6 , implies T k > 0.7, 
consistent with the previous discussion on this parameter 
(Sect. 5.3). 

If more event statistics confirm the large observed op¬ 
tical depths, then the models in Fig. 13a with our best 
COBE-estimate 28° ± 7° for the bar inclination angle are 
inconsistent with the microlensing constraints. This could 
indicate that our models have insufficient mass along the 
line-of-sight in Baade’s Window. To enhance that mass 
without increasing the aq-terminal velocities, an obvious 
possibility is to increase the disc mass outside the bar 
region at the expense of the DH, i.e. approaching a max¬ 
imum disc solution, as suggested by Alcock at al. (1997) 
and by the recent Galactic structure review of Sackett 
(1997). If dark matter in the outer Galaxy is in molecu¬ 
lar form (Pfenniger et al. 1994) and thus concentrated in 
the plane, its contribution to the squared inner circular 


ml2tl600 <?„=19° 



kpc 

Fig. 15. Face-on view of a model with optical depth towards 
Baade’s Window over 3 x IIP 6 . The caption is similar to Fig. 10 
and the model has been symmetriesed 

velocity can become negative , even allowing for an over 
maximum disc solution. To ensure a low central surface 
density, such a heavy disc would then have to be partially 
hollow relative to its exponential extrapolation inside the 
bulge. 

Microlensing optical depths towards Baade’s Window 
also depend on the disc scale height. At fixed surface den¬ 
sity and distance z\ above the plane, one can show that 
the (exponential) scale height which maximises the volume 
density is h z = Z\ if the vertical mass distribution is ex¬ 
ponential. Furthermore, the lenses which most contribute 
to the optical depth are those lying half way between the 
source and the observer and, for reasonable disc parame¬ 
ters, the mass density along Baade’s Window is about con¬ 
stant (e.g. for Hr = 3.5 kpc and h z {R) = const. = 240 pc). 
Thus for bulge sources, the maximum disc contribution in 
this window should come from stars at R = 4 — 5 kpc, 
which are about 270 pc away from the plane. Our mod¬ 
els indeed produce larger values of To on the average 
when h z (R 0 /2) « 300 pc, as shown in Fig. 14. In partic¬ 
ular, models with h z <) 250 pc in the inner part all have 
To < 1.5 x 10~ 6 , hard to conceal with observations. From 
this we infer that the mass distribution of the Galactic 
disc probably has a constant scale height between ~ R 0 /2 
and R 0 , contrary to the constant mass-to-light ratio in¬ 
terpretations of near-IR surface photometry data (Kent 
et al. 1991; Binney et al. 1996). 

The effect of asymmetries like spiral arms may also af¬ 
fect the optical depths towards the bulge, as illustrated 
in Fig. 13a. Model ml2tl600 has one of the largest To to¬ 
wards Baade’s Window and indeed exhibits a strong spiral 
structure (see Fig. 15). Another example, with larger <p a 
and smaller i?L, is m09tl600 (see Fig. 10). The orientation 
of the Galactic bar is such that the line-of-sight through 
Baade’s Window crosses the spiral arm starting at the 
near end of the bar almost tangentially and at a distance 
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Fig. 16a-d. Galactocentric radial kinematical properties of the model rn08t3200. In all plots, the dashed, dotted and full lines 
stand resp. for disc, NS and disc+NS particles, a Radial velocity dispersion and mean velocity within |6| < 5° based on particles 
interior to the Solar circle. The points with error bars are derived from the te Lintel Hekkert et al. (1991) catalogue of double 
peaked OH/IR stars, b Projected radial velocity dispersion along the bar minor axis (1 = 0). The data are from: (1) Sharpies et 
al. 1990; (2) Terndrup published by Rich 1996; (3) Tyson & Rich 1991; (4) Rodgers 1977. c Projected radial velocity moments 
along an axis through the Galactic centre and inclined by 55° relative to the minor axis, and the corresponding Blum et al. (1995) 
M giant observations. The angle a is measured with respect to l = b = 0 and has the same sign as l. d Radial velocity dispersion 
in Baade’s Window as a function of the distance d from the observer, with the K giant data of Lewis & Freeman (1989). The 
lower plot also represents the square root of the relative mass distribution per solid angle along the line-of-sight (arbitrary units) 


where microlensing should be very efficient. The gain in 
optical depth is of order 0.5 x 10 -6 . 

5.5. Inner kinematics 

It is beyond the scope of this paper to study the de¬ 
tailed kinematics of all the dynamical models. Instead, we 
present here the case of a single model, m08t3200, which is 
in fair agreement with almost all explored constraints. In 
Figs. 16a-d, several kinematical predictions of this model 
are compared to stellar observations. 

The agreement between the longitudinal kinematics of 
the model disc and that of the OH/IR stars (Fig. 16a) 
is very convincing. The model particles outside the So¬ 
lar circle were not taken into account as the data con¬ 
tain only very few OH/IR stars at \l\ > 90°. No ve¬ 
locity moment has been derived from the data within 
\l\ < 10° because in this region the radial velocity disper¬ 
sion depends non-negligibly on the galactic latitude and 
the te Lintel Hekkert et al. (1991) OH/IR stars severely 
suffer from undersampling at b <, 2°. 


The minor axis velocity dispersion of the model 
(Fig. 16b) is different for the NS and disc components. 
The observed kinematics are best traced by the NS com¬ 
ponent near the centre and by the disc component at large 
Galactic latitude. In the Blum et al. (1995) fields, the kine¬ 
matics of these components are more similar (Fig. 16c). 

The velocity dispersion in Baade’s Window as a func¬ 
tion of the distance from the observer is close to that ob¬ 
served for K giants (Fig. 16cl), except near the tangent 
point. However, the projected velocity dispersion of these 
stars is obviously less than the 113 km/s of the M gi¬ 
ants which was used to scale the model. Furthermore, the 
distance distribution of the foreground model particles is 
roughly proportional to d 2 , whereas the distribution of the 
K giants seems nearly linear with d (Sadler et al. 1996), 
suggesting that the latter may be biased towards the near 
stars. This difference between model and observations is 
hardly due to a variable Galactic disc scale height. One can 
indeed show that, within 3 — 4 kpc from the observer, the 
simulated line-of-sight mass distribution towards Baade’s 
Window in a realistic analytical double exponential disc is 
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Table 5. Comparison of model m08t3200 versus K giant pro¬ 
jected radial kinematics in off axis bulge fields. Boldfaced, italic 
and roman velocities refer resp. to disc, NS and disc+NS par¬ 
ticles. The mean velocities are Galactocentric. The references 
for the observations are Harding 1996, Minniti 1996a, Minniti 
1996b and Minniti et al. 1992 


in 

b[°] 

[Fe/H] 

Observations 

(J r V r 


Model 

(T r V-p 

-10.0 

-10.0 

> -1 

67 

± 

6 

-82 ± 

8 

72 

-81 



< -1 

107 

± 

6 

-37 ± 

8 

120 

-41 



all 


- 


- 


no 

-71 

9.9 

-7.6 

> -1 

70 

± 

7 

56 ± 

10 

70 

92 



< -1 

91 

± 

13 

-18 ± 

18 

117 

40 



all 


- 


- 


99 

66 

8.0 

7.0 

> -1 

72 

± 

4 

66 ± 

5 

77 

68 



< -1 

109 

± 

10 

—7 ± 

14 

121 

36 



all 


- 


- 


100 

54 

12.0 

3.0 

- 


- 


- 


75 

111 



- 


- 


- 


124 

45 



all 

68 

± 

6 

77 ± 

9 

93 

96 


very similar for a radially constant or a linearly increasing 
h z if the face-on surface density profile is kept the same. 

Table 5 reviews some off-axis K giant observations and 
give the corresponding model predictions. The velocity 
moments of the model disc component resemble those of 
the giants with [Fe/H]> — 1. The l — b proper motion dis¬ 
persions of the model in Baade’s Window are (<r w , a llb ) = 
(3.15,2.38) m"/yr for the disc, (2.96,2.81) m"/yr for the 
NS and (3.08, 2.57) m"/yr for both visible components 
together, whereas the observed values for K giants are 
(3.2 ± 0.1,2.8 ± 0.1) m"/yr (Spaenhauer et al. 1992). 

6. Conclusion 

We have built many self-consistent 3D dynamical barred 
models of the Milky Way extending beyond R = R 0 by 
A-body integration of various bar unstable axisymmetric 
models. The models, extracted from the simulations at a 
frequency of 200 Myr, include 3 components: a nucleus- 
spheroid standing for the Galactic inner bulge and stellar 
halo, a disc mainly representing the Galactic old disc and 
a non-dissipative dark halo. The comparison of the mod¬ 
els with observational constraints leads to the following 
considerations. 

1) The spatial location of the observer in each model is 
constrained by the COBE/DIRBE A-band map corrected 
for extinction by dust, assuming a constant mass-to-light 
ratio Tk for the luminous mass components. The results 
for the best matching models, with mean quadratic resid¬ 
uals between model and data fluxes down to 0.3%, suggest 
that the angle between the l = b = 0 line and the major 
axis of the Galactic bar is 28° ± 7°. 

2) Scaling the models such as the distance of the ob¬ 
server to the centre is R 0 = 8 kpc and the projected radial 
velocity dispersion towards Baade’s Window 113 km/s, as 


observed for M giants, absolute model properties are de¬ 
rived. A dozen of models reproduce fairly well both the 
COBE-data and observations in the Solar Neighbourhood, 
although with a rather low radial versus azimuthal veloc¬ 
ity dispersion anisotropy of the spheroid components. 

3) The bars in these models have a face-on axis ratio 
b/a = 0.5±0.1 and a pattern speed flp = 50±5 km/s/kpc, 
placing the corotation at 4.3 ± 0.5 kpc. Models with a 
disc mass fraction below 0.45 within 3 kpc from the cen¬ 
tre produce envelops of non self-intersecting aq-orbits in 
the l—V diagram which better agree with the observed HI 
terminal velocities. 

4) The microlensing optical depths of the models to¬ 
wards the bulge strongly depends on the bar inclination 
angle <p 0 , increasing from tq ~ 10 -6 for ip 0 = 40° to 
To ~ 2.5 x 10 -6 for ip = 15° towards Baade’s Window, 
whereas observations rather support values over 3 x 10 -6 . 
We find that a spiral arm starting at the near end of the 
bar can increase the optical depths by 0.5 x 10 -6 and thus 
reduce the gap between model and observed values. 

5) All models with a disc scale height h z < 250 pc half 

way between the observer and the Galactic centre have 
To < 1.5 x 10 , arguing against an inwards decreasing 

disc scale height. This result is also in agreement with the 
constant h z observed in late-type spirals. 

6) The models predict a mass-to -K luminosity ratio 
T k = 0.6—0.8 in Solar units. Values near the upper limit 
are consistent with the mass-to-light ratio estimated for 
M100 (NGC4321), a galaxy similar to the Milky Way, and 
favour microlensing optical depths closer to the observed 
values. There is indeed an obvious but tight correlation 
between To and T k, according to which To <1 2 x 10 -6 
implies T k > 0.7. 

7) Most models predict too low surface brightness rel¬ 
ative to the COBE-data in the region |/| ^ 15° dominated 
by the disc. Beside an improbable lower disc scale height in 
the inner Galaxy, or a variable mass-to-light ratio, this dis¬ 
crepancy could indicate that the Galactic disc outside the 
bar region is more massive than assumed in the models, 
favouring large microlensing optical depths and arguing 
for a maximum disc Milky Way. 

8) The disc radial kinematics of the models resembles 
the observed kinematics of K giants with [Fe/H] > — 1 in 
the outer bulge. 

As reasonable models regarding most observational 
constraints, we would recommend the models m08t3200 
and m04t3000. 
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Appendix A: a non-Gaussian bounded 3D distri¬ 
bution 

A convenient 3D distribution with non-zero probability 
over a finite volume, i.e. avoiding the tails of the multi¬ 
normal distribution, is given by: 

B 3 (Z,v, 0 = ■ 


(1+0 (1-0 (!-T 


r) M 2 (I - 


VCD— 1 


if 0 + ? 7 2 + C 2 < 1 


(Al) 



0 


otherwise, 


where B is the Beta function, and k. A, p and ui are four 
parameters. In the reduced variables £, 77 and £, this dis¬ 
tribution is bounded by a sphere of radius 1 on which, as 
long asK, A>2,/r>3/2 and u> > 1, it continuously van¬ 
ishes. If k = A = 2, n = 3/2 and u> = 1, the distribution 
is homogeneous inside the boundary sphere and smaller 
values of these parameters produce singularities on it. 

The first and second moments are: 


v = C = (£ - = (Z - 0C = vC = 0 , 
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Depending on k and A, the distribution is skewed in £, the 
maximum of probability lying at: 


Fig. Al. Distribution of the variables f and 7 resulting from 
Eq. (A 8 ), with k — 5, A = 3 and 77 = 4 
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To generate random numbers distributed according to B 3 , 
one may use the property that: 

k = 2 ^ + 1 )’ 

(A13) 

t - l ( V + 11 , 

v 2 \/i-e 

(A14) 

k ~dv^e^ +1) 

(A15) 


follow Beta distributions with parameters k. A for t%, 7,7 
for t v , and u>,u> for fy. 
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Appendix B: the TZ 2 residual 


Figure Al displays an example of the 2D distribution B 2 
obtained after integrating B 3 over Q: 

= J B 3 (£,r],C)d( 


If the model fluxes are decomposed into F,; = F., \ + Fj j2 , 
where F it i is the exact part and Fj j2 is the statistical error 
due to the finite number of particles, the y 2 in Eq. (16) 
expands into: 


(1 + 0 *— 3 / 2 ( 1 - 0 A “ 3 /2 ( 1 - T ^ 2 ) • 


(AS) 


The inversion of Eqs. (A3)-(A6) provides the param¬ 
eters of the ^ 3 -distribution as a function of the aimed 
moments: 


iVpi* 
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(Fi, 1 - F °) 2 


N pix p 2 

h a i h a i 


(Bl) 


i(i+D(i-f -4) 
2 4 

i(i-0(i-f -4) 
9 ^2 


, . . The expected value of the second term in the right hand 
' side, when resulting from a fit, is the number of degree 
of freedom v (if the errors were Gaussian, then this term 
(A10) wou ld follow the standard y 2 statistics), whereas the last 
term is about zero since by definition the F (j2 ! s must van- 
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ish on the average. With these simplifications, the best 
estimate of the function: 


n 2 = 


N pix 

E 


o*/F° 


N plx 

E 


«i/f ° 2 


(B2) 


i.e. the quadratic relative residuals between the F lt i's and 
the F° ’s averaged over pixels and weighted by the inverse 
of the relative variance, indeed reduces to Eq. (15). 


Appendix C: the variance of the model fluxes 

The model flux in a given pixel is of the form: 

N 

F{D) = Y J f{D k ), (Cl) 

fc= 1 

where N is the total number of particles in the pixel, the 
D k s are the distances of the particles relative to the ob¬ 
server and D = (Hi,...,Djv). Noting resp. P(N) and 
p{D) the distributions of N and of the distances, the two 
first moments of F are: 

°° nO O POO 

F^=y^P(N)-l dD lP (Di)... dD N p(D N )F(D) n 

Jo 

n = 1 

(C2) 

n = 2, 


(C3) 

n = 1,2. (C4) 

The result for n = 2 in Eq. (C2) requires the assump¬ 
tion that the first and second moments of P are identical, 
which is true for Poisson statistics. Hence the variance of 
F becomes: 
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